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A global groundwater flow model has been constructed for Mars to investigate hydrologic 
response under a variety of scenarios, improving and extending earlier simple cross-sectional models [1]. 
The model is capable of treating both steady-state and transient flow as well as permeability that is 
anisotropic in the horizontal dimensions. A single near-surface confining layer may be included 
(representing in these simulations a coherent permafrost layer). Furthermore, in unconfined flow 
locations of complete saturation and seepage are determined. The flow model assumes that groundwater 
gradients are sufficiently low that DuPuit conditions are satisfied and the flow component perpendicular 
to the ground surface is negligible. In spherical coordinates the governing equation is: 
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where X is longitude, d> is latitude, h is the local hydraulic head, R is the planetary radius (assumed 

constant), S is the specific storativity (confined portions of the aquifer) or specific yield (unconfined 

areas), Q is distributed vertical recharge, and are the three components of the symmetric 

hydraulic conductivity tensor, and he is an effective aquifer thickness which, in the case of an aquifer with 
K that does not change vertically (termed the uniform case) is equal to (h u - h\), where h u and hi are the 
elevations of the top and base of the aquifer, respectively. In confined situations h e equals (h c - hi), 

where h c is the base of the confining layer. In some simulations K was assumed to decrease 

exponentially with depth from the surface at elevation h s (the exponential case), giving an effective 
thickness of 
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where a is a decay constant (in these simulations the aquifer is assumed to be unbounded at depth). In 
these simulations storativity was assumed also to decrease exponentially with depth with the same decay 
constant. In unconfined regions h u equals the water table elevation, but in confined circumstances h u is 
the base of the confining layer, h c . 

Surface elevations were taken from the new 1:1 5M Mars topographic map. Two sets of elevations 
were recorded for each grid-point location, the first (the nominal set) being the local elevation, which, 
however, was adjusted upwards to the level of the surrounding uplands if the gridpoint occurred within a 
crater or in an erosional channel. However, locations in craters large enough to encompass more than one 
grid points (e.g., Hellas and Argyre) were represented by the local elevation. The second set of 
elevations (the fluvial set) was taken in the same manner as the nominal set except that the elevation was 
adjusted downwards to the elevation of the bottom of any nearby erosional channel. Thus the nominal set 
is assumed to represent the pre-channeling landscape and the fluvial set to represent the present fluvial 
base level. Examination of both outflow and the larger cratered-terrain valley networks reveals that most 
flowed down essentially the present regional topographic gradients. The few local disparities between 
inferred channel flow directions and topographic maps (e.g. Nirgal, Dueteronilus, Mangala) are probably 
as likely to be due to uncertainties in topographic mapping as to regional tectonism. The major exception 
is the Tharsis volcanic construct, largely postdating the development of valley networks. However, flow 
directions of most outflow channels associated with the margins of Tharsis appear to be consistent with 
present topographic gradients, suggesting that most of the marginal updoming and normal faulting in 
Valles Marinaris predated the outflow channels. In addition, the direction and degree of development of 
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the valley networks suggests that the north-south topographic discontinuity (highlands-lowlands transition) 
had essentially its present configuration prior to the development of the late Noachian valley networks. 
Thus the present topography reasonably can be considered to be indicative of conditions during channel 
development. 

In simulations with a finite aquifer thickness, the local base of the aquifer, h\, was set to a given 
elevation below the local fluvial elevations. In the present simulations the aquifer thickness below the 
fluvial elevations was assumed to be areally uniform. 

The flow equations were solved using a finite difference method employing 10-degree spacing of 
latitude and longitude. A successive over-relaxation method (SOR) was used for steady-state solutions 
and a Crank-Nicholson method with iterative solution at each time step was used for transient simulations. 
The correctness of the solution method was checked in part by confirming mass conservation. In 
addition, an independent, finite element flow model was also constructed for steady-state conditions and 
compared to the finite difference model. The finite element model represented the planet by triangular 
planar elements matched at their edges and the flow equations were cast in cartesian form. The two 
models gave essentially identical solutions. 

The initial guess for steady state flow simulations was a level water table with an elevation equal 
to the lowest surface elevation on the planet (the bottom of Hellas). If the base of the aquifer was above 
this level then the water table was set to the base of the aquifer. Steady state simulations assumed a 
constant recharge of the aquifer from the surface. In most cases an areally uniform recharge rate was 
assumed, but some simulations were conducted that assumed recharge to be a function of latitude and/or 
elevation: 


Q = Qo [CosO ] {1 + /3 A s ) 

where £>o is a nominal recharge rate, and /? is an input parameter. The bracketed latitude correction is 
optional. Recharge to fixed-head and confined portions of the aquifer was set to zero. The only location 
on the planet that was specified as fixed-head at the beginning of the iterative solution was the bottom of 
Hellas. Locations of saturation and, therefore, seepage were determined during the SOR iterative solution 
by identifying during each iteration locations where the predicted head, h, was greater than the surface 
elevation, h s . Such locations were converted to fixed-head for succeeding iterations until the steady-state 
solution was achieved. A similar procedure was used to identify confined portions of the aquifer 
wherever the predicted head was greater that the bottom of the confining layer, h c . 

Transient flow simulations started from steady state conditions, and involved either of two 
scenarios: 1) Draining : The aquifer system is initially equilibrated with a specified recharge rate and 
then drained with no further recharge, and 2) Filling: The aquifer system initially has no available water 
(i.e., it is set to the initial conditions for the steady state iterations as described above) and the aquifer 
then fills towards steady state with a specified constant recharge. 

Steady-state flow simulations. Steady-state flow simulations have been conducted to investigate 
effects of parameter variation, particularly the ratio of recharge rate, Q y to intrinsic permeability, k . 
Hydraulic conductivity has been assumed to be isotropic and uniform in these simulations. Another set of 
runs examines the effects of a confining permafrost layer extending from the poles to a variety of 
latitudes. Finally some simulations have been conducted to examine the effects of spatially variable 
recharge and anisotropic permeability. 

Transient flow simulations. A two-dimensional flow model had previously been used to calculate 
the length of time for martian aquifers to drain following cessation of recharge []]. These figures are 
updated here using the global flow model. Results from a number of simulations run with different sets of 
model parameters were analyzed by multiple regression to determine estimating equations for the length 
of time for specified percentages of filling, tf> or draining, t& of the aquifer. The estimating equation for 
draining flows is given by 
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where C<j is a coefficient (see table below), v is the water viscosity, tj is the porosity, k is the intrinsic 
permeability (note that hydraulic conductivity, K equals kg/v), g is the gravitational constant, G is the 
average hydraulic gradient at the start of aquifer draining, and V t is the ratio of the available water 
volume (water above the lowest discharge point in Hellas) to the total water volume (including water 
below the lowest discharge point) at the start of aquifer draining. The exponents take the values of 6=-1.3 
and y=l for the uniform case and &=-1.0 and )t-1.4 for the exponential case. Similarly, for the case of 
filling flows: 
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where D is the depth from the surface to half-value of the permeability (0.693/a) and the exponents tak*. 
the values of o^-0.8, e=0.2 and y=0 for the uniform case and a^-0.8, e=0.2 and y=\.3 for the exponential 
case. The times to filling or draining can be specified in two ways: 1) the time required to drain or fill a 
specified percentage of the total available water volume (constants Q v and Cfv); or 2) the time required 
for the total discharge of water to decrease or increase by a specified percentage (constants Cd q and Cf q ). 
The values of the coefficients are given below: 
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In using this table Q is expressed in cm/yr, D in km, k in darcies, g in cm/sec^, v in cm^/sec, and times, t, 
in 10 6 yr. 


The times for aquifer draining are generally comparable to those reported for the two-dimensional 
simulations [1], For example, for the uniform permeability case with tj=0.2, k=l, G- 0.0017 and VV=0.67 
the 2-D simulations with a representative aquifer length of 3000 km indicate the time to 75% draining as 
5.6*106 yr, and the present simulations indicate 12X10 6 yr for the uniform aquifer and 48*10® yr for the 
exponential aquifer. However, the draining times increase much more strongly with percent draining than 
in the case of the 2-D simulations because the length-scale increases as draining progresses and the water 
must flow towards more distant exit points. 

Times for aquifer recharge are essentially inversely proportional to recharge rate. Furthermore, 
recharge generally occurs much more rapidly than draining for comparable regolith permeability. For 
example, recharge at a rate of 1 cm/yr will fill the regolith to 75% of its steady-state capacity in only 
150,000 years (assuming 7j=0.2 and k- 1) compared to tens of millions of years for a comparable 
percentage of draining. In addition, recharge rates are only weakly dependent upon aquifer permeability. 
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